Model of Magnetically Shielded Ferrite-Cored Eddy Current Sensor

Computationally fast electromagnetic models of eddy current sensors are required in model-based measurements, machine interpretation approaches or in the sensor design phase. If a sensor geometry allows it, the analytical approach to the modeling has significant advantages in comparison to numerical methods, most notably less demanding implementation and faster computation. In this paper, we studied an eddy current sensor consisting of a transmitter coil with a finitely long I ferrite core, which was screened with a finitely thick magnetic shield. The sensor was placed above a conductive and magnetic half-layer. We used vector magnetic potential formulation of the problem with a truncated region eigenfunction expansion, and obtained expressions for the transmitter coil impedance and magnetic potential in all subdomains. The modeling results are in excellent agreement with the results using the finite element method. The model was also compared with the impedance measurement in the frequency range from 5 kHz to 100 kHz and the agreement is within 3% for the resistance change due to the presence of the half-layer and 1% for the inductance change. The presented model can be used for measurement of properties of metallic objects, sensor lift-off or nonconductive coating thickness.


Introduction
Eddy current sensors (probes) respond to changes in the material properties (electrical conductivity, magnetic permeability), geometry or position of the nearby electrically conductive, usually metallic, objects [1,2]. Their applications include nondestructive testing for defects and material changes in objects, such as heat exchangers in nuclear plants, oil well casings, pipelines or aircraft hulls [3][4][5]; measurement of the distance to such objects or their dimensions [6,7]; and material characterization [8][9][10]. Eddy current sensors offer detection of very small surface defects, measurement of dimensions or distances in the micrometer range, insensitivity to nonconductive dirt, such as oil, water or dust, and robustness to the environmental factors such as pressure, temperature or stress. However, absolute measurement using eddy current sensors can be hard to achieve leading to the need for calibration, and human or machine analysis and interpretation of the sensor signals [11][12][13].
Traditionally, an eddy current sensor consists of a transmitter coil that generates a magnetic field and induces the eddy currents in the conductive material [1]. Receiver coils pick up the resultant magnetic field due to the defect and changes in the dimensions and material properties. Receiver coils can be replaced by Hall, magnetoresistive or even superconducting quantum interference devices (SQUID) sensors [14]. However, the overall characteristics of the eddy current sensors are primarily determined by the excitation, geometry and construction of the transmitter coil.
Fast electromagnetic models of transmitter coils are required in model-based measurements, machine interpretation approach or in parametric and sensitivity analyses in the sensor design phase. When it comes to more complicated sensor geometries, analytical electromagnetic models are not as versatile as numerical ones, e.g., finite-element or boundary-element models. However, if a real sensor and object under test can be approximated by a geometry that is simple enough to be dealt with analytically, then one can expect less demanding numerical implementation and higher computational speed in comparison to fully numerical approaches. More complicated geometries can be treated using hybrid approaches in which the analytical models are used for calculation of the incidence fields and numerical methods for calculation of the fields scattered from the small changes in the tested object [15].
Typically with eddy current sensors, the transmitter coil has a ferrite core to increase the inductance, sensitivity and signal-to-noise ratio or to achieve the desired magnetic field distribution. Since the seminal work by Dodd and Deeds [16], who studied air-cored transmitter coils in planar and cylindrical geometries (layers, rods, tubes), a number of more complicated geometries have been analyzed using the truncated region eigenfunction expansion approach (TREE) [17]. This includes ferrite-cored transmitter coils, which previously required numerical or empirical approaches. In a paper that laid out the modeling path followed by other authors, Theodoulidis analyzed the coil with the cylindrical finite-length ferrite core that was positioned above a half space made of two nonmagnetic and conductive layers [18]. This analysis was extended to a half space made of multiple magnetic and conductive layers by Lu et al. in [19]. Bayani et al., in [20], developed the model of a cup-cored coil above a nonmagnetic conductive double layered half space. In all three studies, the models were verified experimentally and the relative error between the theoretical results and measurements were on average between 0.5% and 2.5%. Sakkaki and Bayani, in [21], modeled an E-cored coil above a conductive double layer, and they numerically verified the model by comparing it to a finite element model using Comsol. Tytko and Dziczkowski, in [22][23][24], modeled the I-cored and E-cored coils in the presence of the layered medium with a hole, and they verified their work using a Comsol-based finite element model. Zhu et al.,in [25], analyzed an I-cored coil screened with a thin sheet of material with infinite permeability.
In this paper, we studied a coil with a finitely long I core and screened with a finitely thick magnetic shield placed above a conductive and magnetic half-layer using the TREE approach. Figures 1 and 2 depict the geometry in more detail. We made no assumptions regarding the shield thickness and the permeability of the shield and the core. The open upper side of the shield allows for the coil wires and connector towards a measurement instrument. This sensor geometry is typical and practically important, but to the best of our knowledge it was not modeled analytically before. The contributions of this paper are the analytical model of the sensor impedance above a plate made of conductive and magnetic material in Section 2 and accompanied Appendix A, analysis of numerical implementation issues in Section 3, comparison with a finite element model in Section 4, experimental verification in Section 5 and sensitivity analysis with respect to the sensor dimension in Section 6. The immediate application of these contributions could allow model-based measurement of material properties and lift-off (e.g., nonconductive coating thickness) with no or less calibration points (measurement above reference standards with known material and lift-off) compared to the existing measurement procedures [26]. The model can also be used in the sensor design phase where it is often necessary to adjust the sensor properties to the specific application or account for sensor production tolerances. Furthermore, we believe that the modeling approach and discussion on the implementation issues can be useful in modeling other similar induction sensors with ferrites.  A ferrite-cored coil screened with a ferrite shield above magnetic and conductive half space. A filamentary coil is also shown. For better clarity, the subregions of a region are indicated by colors with the same hue value but different lightness. In the case of the filamentary coil, region C is not present and the boundary between regions 3 and 4 is at z = z 0 .

Analysis
The geometry of the problem is shown in Figures 1 and 2. The coil with a rectangular cross section is wound around the ferrite core with relative magnetic permeability µ C . The coil and its core are screened with the ferrite shield with relative magnetic permeability µ S . The core and the shield are vertically aligned at the plane z = 0. The core and the shield are assumed to be nonconductive. The sensor is located over the magnetic and conductive half space (plate) with electrical conductivity σ and relative magnetic permeability µ r .
Before we derive the vector magnetic potential of the coil, we will first analyze a filamentary sinusoidal current I exp(jωt)δ(r − r 0 )δ(z − z 0 )φ, where δ is the Dirac delta function andφ is the azimuthal unit vector for the cylindrical coordinate system. The problem is axially symmetric, so there is only ϕ component of the vector magnetic potential, i.e., A(r, z) = A(r, z)φ.
The problem domain is truncated at far enough boundary r = R, where we will impose the Dirichlet boundary condition A(r = R, z) = 0. The consequence of the truncation is that the final expressions for the potential are in series rather than integral form. The problem domain is divided into horizontal regions indicated generally by index n in the expressions, or specifically by a number or the letter C. If there are multiple materials or the coil in a region, the region is divided into subregions that are indicated generally by index m or specifically by a number following the region designation. For example, A 1 is the potential in region 1 (no subregions) and A 23 in region 2, subregion 3.
Following the separation of variables, the general form of the magnetic potential is U n,i J 1 α n,i r + V n,i Y 1 α n,i r C n,i exp −β n,i z + D n,i exp β n,i z , for n ∈ {1, 5, 6}, where J 1 and Y 1 are first-order Bessel functions of the first and second kind, respectively, α n,i are eigenvalues arising from the Dirichlet boundary condition at r = R for each region, Section 2.2, and If a region is homogeneous with relative permeability µ n and conductivity σ n , i.e., it has no multiple subregions, the form of the potential is given by (1). For regions with multiple subregions, the general form of the potential in subregion m of region n with relative permeability µ nm is given by (2). It is important for further steps to note that there are no conductive materials in the regions with multiple subregions (only air or ferrites) and that, consequently, both r-related and z-related parts have the same eigenvalues α n,i . This means that z-related parts of the subregions in a given region are the same in case of (2). The homogeneous regions n ∈ {1, 5, 6} have their corresponding coefficients V i = 0 because of the divergence of Y 1 at r = 0, so, consequently, one can freely set U i = 1. Coefficients U nm,i and V nm,i in (2) for regions n ∈ {2, 3, 4, C} are determined from the interface conditions along the radial boundaries between the subregions, as described in Section 2.1. The interface conditions in this case are continuity of B r and H z , i.e., component of B inr direction, and component of H inẑ direction.
The unknown coefficients C n,i and D n,i in (1) and (2) are determined from the boundary and interface conditions along the horizontal boundaries between the regions, Section 2.4. The potential must remain finite as z → ±∞, so D 1,i = 0 and C 6,i = 0. The interface conditions for the horizontal boundaries are continuity of B z and H r .

Continuity of B r and H z at Radial Boundaries between Subregions
This condition has to be satisfied for regions 2, 3 and 4. We will analyze a general case of a region n with subregions m = 1, . . . , M and radial interfaces between these subregions at r n1 , . . . , r n(M−1) . In that case, the radial parts under the sum in (2) for each subregion are R n2 (α n,i r) = U n2,i J 1 (α n,i r) + V n2,i Y 1 (α n,i r), r n1 < r ≤ r n2 , (5) . . . R nm (α n,i r) = U nm,i J 1 (α n,i r) + V nm,i Y 1 (α n,i r), r n(m−1) < r ≤ r nm , . . . R nM (α n,i r) = U nM,i J 1 (α n,i r) + V nM,i Y 1 (α n,i r), r n(M−1) < r ≤ R, r nM = R.
For the two neighboring subregions m and m + 1, the radial interface of continuity of B r and H z at r = r nm is equivalent to Since all subregions in region n have the same z-related part, the radial interface conditions can be satisfied by adjusting the coefficients of the r-related part only. The notation will be simpler if we define matrix R as It can be shown using (2)-(10) that the coefficients in the r-related part of subregion m + 1 are connected to the coefficients of subregion m as Using recurrence relation (11), and since V n1,i = 0 and U n1,i can be freely set to 1, we can determine the r-related coefficients for all subregions in a given region: for m = 3, · · · , M, r n(m−1) < r ≤ r nm , r nM = R.
Because coefficients U nm,i and V nm,i do not depend on spatial coordinates r and z, the radial parts R nm (α n,i r) given in (4)-(7) are linear combinations of Bessel functions J 1 and Y 1 with respect to r. This will be important later on when we make use of the orthogonality property of the Bessel functions. In order to facilitate this, we will add an index to R nm to indicate the order of the involved Bessel functions: where coefficients U nm,i and V nm,i do not depend on ν. It follows from (8) and (9) that

Eigenvalues
The radial parts of all regions and subregions are defined in the preceding section. The first step in their computation is to find the eigenvalues that satisfy the Dirichlet condition at the truncated boundary of the problem: The eigenvalues α n,i are the real and positive roots of (18) and (19). Homogeneous regions 1, 5 and 6 have identical radial parts and, hence, the same eigenvalues. Similarly, the regions 3, 4 and C, with multiple subregions, have the same radial parts and the same eigenvalues. Region 2 has its own eigenvalues. In order to simplify the notation, which was so far useful for the general case, we will designate the elements of the three sets of eigenvalues as α i , p i and q i : and for the z-related part in region 6

Final Expressions for Potential in All Regions in Matrix Form
Infinite series in (1) and (2) have to be truncated to N terms in the computation. The potential is In the above expressions, we divided with the corresponding eigenvalues in order to make later expressions somewhat simpler, which is allowed because it will be canceled out by the final form of C n,i and D n,i [17]. It is convenient to write (24)-(30) in the matrix form. To do that, we will assume that the functions J ν (), R ν,nm (), exp(), as well as the integration and differentiation are applied element wise on a vector or matrix. We will use boldface type to represent matrices (e.g., α or exp(qz)), overlined italic type to represent row vectors (e.g., α or J 1 (αr)) and underlined italic type to represent column vectors (e.g., C 1 ). The dimensions of the matrices are N × N, row vectors 1 × N and column vectors N × 1. This notation eases implementation in a programming language such as Matlab or Julia. Using this notation, the matrix form of (24)-(30) is In the above expressions, matrices α, p and q are diagonal with elements of α, p and q on their main diagonal, respectively.

Continuity of B z and H r for Horizontal Boundaries between Regions
The continuity of B z and H r have to be satisfied along the horizontal interface (z = const.) between two regions. For two neighboring regions n and n + 1, the horizontal interface condition at z = z n is equivalent to Relative permeability µ n in (38) is the piecewise constant function of r for region n if it has multiple subregions, i.e., µ n = µ n1 for 0 ≤ r ≤ r n1 ,. . . , µ n = µ nM for r n(M−1) < r ≤ R.
The interface relations that have to be satisfied for all five horizontal boundaries (z = L S , z = L C , z = z 0 , z = 0 and z = −h) are given in Appendix A for the sake of clarity. In order for these relations, e.g., (A3), to be satisfied for the entire interval 0 ≤ r ≤ R, the procedure followed in Appendix A is to expand the r-dependent parts using a basis formed of continuous piecewise linear combinations of Bessel functions of the first and second kind. The choice of the basis is somewhat arbitrary but it is better, from the standpoint of numerical implementation, to choose these in such a way that as many of the resulting matrices are identical or diagonal. Term-by-term comparison results in a matrix equation for each interface relation involving corresponding vectors C and D. This procedure has to be repeated twice (continuity of B z and H r ) for each of the 5 horizontal interfaces giving rise to 10 algebraic equations with matrix coefficients and 10 unknown vectors C and D.

Calculation of Coefficients C and D
For the filamentary coil at (r 0 , z 0 ): For the coil with the rectangular cross section: where int(x 1 , Finally, we can write the coefficients C and D sorted in order of their computational dependencies:

Coil with Rectangular Cross Section
The magnetic potential of a coil with rectangular cross section and N TX turns, carrying sinusoidal current of amplitude I TX , is calculated by the superposition of a number of filamentary coils. Let a coil with infinitesimal cross section dr 0 dz 0 carrying current i TX dr 0 dz 0 approach a filamentary current I at (r 0 , z 0 ), where the current density i TX is In that case, the total magnetic potential for region n is where A n (r, z, r 0 , z 0 ) is given by (31)-(36) depending on the region. Integration in (66) comes down to integration of X pos X and X neg X in (48)-(50), and it results in (51)-(53). The total potential A tot n is given then by the same expressions as A n in (31)-(36), except one uses (51)-(53) instead of (48)-(50) in (55), and (60) and (61) for D 5 , C 3 and D 3 , respectively. One should note that superposition for region C (z T1 ≤ z ≤ z T1 ), which is present only in the case of the rectangular cross section coil, requires integration of A 3 (r, z, r 0 , z 0 ) for z T1 ≤ z 0 ≤ z and A 4 (r, z, r 0 , z 0 ) for z ≤ z 0 ≤ z T2 , i.e., It can be shown that (67) results in where X is given by (51). The impedance of the coil, which is the most important output of the model from the measurement point of view, is calculated by the integration of the voltage induced in a single loop over the coil cross section: where A C2 is the potential in region C, subregion 2. Finally, the impedance of the coil is The voltage induced in a receiver coil can be calculated using an approach analogous to (69) with the choice of the potential corresponding to the location of the receiver coil.

Numerical Implementation
The model can be implemented in a programming language suitable for numerical analysis, such as Matlab or Julia. We have chosen the former. The major functional parts of the implementation are: 1. Calculation of coefficients for the radial part R nm (α n,i r) given in (4)-(7) using (12)-(14). 2. Root finding of (18) and (19) in order to obtain the required eigenvalues. 3. Integration of int(x 1 , x 2 ) in (54). 4. Calculation of the matrices in Appendix A (E, K, F, W, L, U, V and S) to form the set of 10 algebraic equations for 10 unknown coefficients C and D. 5. Calculation of the coefficients C and D using equations in Section 2.5. 6. Calculation of the coil impedance Z using (70).
While steps 1, 4 and 6 were implemented straightforwardly following the relevant expressions in the paper, the remaining ones deserve some additional comments.
The functions on the left side of Equations (18) and (19) are real and oscillatory. We found their roots by sampling the functions so that intervals containing exactly one zero crossing can be selected and the roots bracketed. As a root finding algorithm, we used Anderson-Björck modification of the classic regula falsi method [27].
In general, the form of the integral in (54) involving linear combination of Bessel functions of the first and second kind, J 1 and Y 1 , is more easily calculated using Struve functions or expansion in Chebyshev polynomials than by direct numerical integration [18]. We used Struve functions H 0 and H 1 calculated via their expansion into the series of Bessel functions for which the relevant identities can be found in [28]: The expressions for C and D in Section 2.5 require some matrices to be inverted. Increasing the length N of the series representation increases the eigenvalues. The factors containing exponential functions with positively signed products of an eigenvalue and a spatial dimension can cause a loss of precision and a poor condition number of the matrices that need to be inverted. Thus, care should be taken to arrange the final expressions in such a way that as much as possible of the matrices are diagonal and that the exponential functions have arguments with negative sign.

Comparison with FEM
In order to verify the model of the sensor impedance, we compared the model predictions with the results of a finite elements model (FEM) study. In both numerical studies we used the sensor with properties as in Table 1. For the FEM study we used free, simple and thoroughly verified 2D solver Finite Element Model Magnetics (FEMM) that was successfully used in a number of studies [29]. The number of the series terms in the presented model of the sensor impedance was N = 140 and the domain boundary was set to R = 10r S2 . The identical geometry in FEMM was meshed with approximately 2.4 million of triangular elements. The computation time for one realization of the problem (i.e., for the selected dimensions, frequency and properties of the core, shield and medium) was around 30 min for FEMM, and less than 0.7 s for the calculation of the model, including the eigenvalues computation and without usage of precalculated values.
First, we calculated the impedance of the coil depending on the permeability of the core and shield for the plate with σ = 5 MS/m and µ r = 50 at 60 kHz, Figure 3. For the sake of simplicity, we assumed that µ C = µ S . There is practically no significant further change in the impedance of the sensor for the permeability above 500, which is a finding observed for other cored coils as well [17]. The relative discrepancy between the model and FEM results are bellow 0.07% for the inductance and 0.23% for the resistance. These discrepancies can be further reduced by increasing the density of the FEM mesh.
Second, we calculated the impedance of the coil depending on the conductivity σ of the plate for two values of its permeability µ r = 1 and µ r = 50 at 60 kHz, Figure 4. The core permeability µ C = 100 and the shield permeability µ S = 50. Typical dependency of the coil impedance on the plate conductivity can be observed, including the maximum point of the resistance [17]. The relative discrepancy between the model and FEM results are bellow 0.07% for the inductance and 0.25% for the resistance.
Finally, we calculated the impedance of the coil depending on the lift-off h of the sensor placed above either a nonmagnetic plate with µ r = 1 and σ = 25 MS/m (properties similar to an aluminum alloy), and a magnetic plate µ r = 50 and σ = 5 MS/m (properties similar to a carbon steel), Figure 5. The core permeability was µ C = 100 and the shield permeability µ S = 50. Again, the relative discrepancy between the model and FEM results are bellow 0.07% for the inductance and 0.25% for the resistance. As expected, the resistance decreases with the lift-off increase for both plates, whereas the inductance in the case of the nonmagnetic plate increases with the lift-off, and decreases in the case of the magnetic plate.

Comparison with Experimental Results
We manufactured the sensor according to the geometry in Figures 1 and 2. The sensor properties are given in Table 2. The impedance of the sensor in the air and above an aluminum plate was measured using the precision LCR meter HP 4284A in the frequency range from 5 kHz to 100 kHz and for lift-off values from 0 µm to 1000 µm in 100 µm steps. The plate was made of aluminium alloy EN AW-6082 with a conductivity of 42.2% IACS or 24.56 MS/m. The impedance of the coil Z is the sum of the coil impedance in air (without the plate) Z 0 = jωL 0 , and the impedance change ∆Z = ∆R + jω∆L due to the proximity of the plate: where the resistance of the coil in air and the coil parasitic capacitance have been compensated for. The inductance of the sensor in the air was calculated 1525.67 µH, while the measured value was from 1524.48 µH to 1526.38 µH depending on the frequency. The impedance change obtained from the model and the experiment are shown in Figures 6 and 7 for ∆R and ∆L, respectively. The model and experiment are in practically acceptable agreement with the relative error of ∆R within ±3% and ∆L within ±1%. These errors are comparable to the ones in other studies on eddy current sensor modeling. In practice, model-based measurements of conductivity and lift-off require calibration procedures with samples of known conductivity and lift-off if one wants to achieve measurement uncertainty better than the reported modeling errors would allow.

Sensitivity to Sensor Dimensions
The data in Table 2 was obtained by averaging and rounding to 0.01 mm the values measured using a Mitutoyo micrometer with resolution of 0.001 mm and accuracy of ±2 µm. We rounded the values because at a higher resolution, the measurement is affected by the applied pressure or measurement position along the measured object, e.g., the ferrite core was not a perfect cylinder at a resolution of 0.001 mm. We used purposefully overestimated measurement uncertainty of ±0.01 mm to evaluate its effect on the sensor impedance. With that aim, we calculated the sensitivity of the sensor impedance to the sensor dimensions. Table 3 shows the relative changes of the inductance in the air L 0 and impedance change ∆R and ∆L (for plate with σ = 5 MS/m and µ r = 50 and lift-off 0 µm) if each of the nine sensor dimensions is independently varied within the interval of ±0.01 mm. The core radius and the shield inner radius uncertainties have the most significant contribution to the impedance uncertainty. This means that the gap between the core and the shield must be precisely controlled during the sensor production. Table 3. Relative uncertainty of sensor impedance in the case of ±0.01 mm uncertainty in sensor dimensions. The plate has σ = 5 MS/m and µ r = 50, lift-off is 0 µm and the frequency is 60 kHz.

Conclusions
The validity of the impedance model is corroborated by a practically negligible discrepancy in comparison with the FEM approach. We compared the model and FEM approach for a range of practical scenarios: changes of the impedance due to the variations in the core and shield permeabilities, dependence of the sensor impedance on the plate conductivity of magnetic and nonmagnetic plates and impedance dependence on the sensor lift-off above magnetic and nonmagnetic plates. The numerical implementation of the model is very fast in comparison to the FEM approach (0.7 s vs. 30 min), and this can be improved further if the intermediate results are reused for the next iteration. In comparison with the impedance measurement, the relative error better than 3% in the frequency range from 5 kHz to 100 kHz confirms the applicability of the model for eddy current measurement of material properties and lift-off. The model can be extended easily to include multi-layered planar structures.

Institutional Review Board Statement: Not applicable.
Informed Consent Statement: Not applicable.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author.

Conflicts of Interest:
The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses or interpretation of data; in the writing of the manuscript or in the decision to publish the results.

Appendix A. Interface Relations for Horizontal Boundaries
The following integrals involving the Bessel functions will be very useful: where X and Ψ are any two cylinder functions.
Because of the piecewise definition of the radial part of the potential, the following equations are obtained from the continuity of B z according to (37): If both sides of (A3) are multiplied from the left by the basis function, and then integrated from 0 to R, we obtain where elements of matrices E = (e ik ) and K = (k ik ) are where p i and α k are elements of the vectors p and α, E is full, and K is diagonal. Continuity of H r according to (38) results in the following: Appendix A.2. Interface z = L C Continuity of B z according to (37) results in the following: The basis function X pr is given by (A4).